$ONTEXT
David Morton
	A capacity expansion planning model under uncertainty. There are two types
	of uncertainty: (i) demand and (ii) generator availability. In the first
	stage, capacity is allocated to two generators. In the second stage, a
	transportation problem is solved to send electricity from the generators to
	three demand sites. If the available capacity is not sufficient to meet the
	demand then we must purchase electricity from an outside source (subcontract)
	to satisfy the demand.
$OFFTEXT
$TITLE capacity expansion planning problem
$OFFUPPER OFFSYMXREF OFFSYMLIST
$INLINECOM { }

OPTIONS ITERLIM = 10000, RESLIM = 10000, LIMROW = 0,
	LIMCOL = 0, SYSOUT = OFF, SOLPRINT = OFF,
	LP = cplex, OPTCR = 0.00001;
SETS
	I generators /I1*I2/
	J demand sites /J1*J3/
	A1 available from generator 1 /F1*F3/
	D1 demand at site 1 /L1*L3/;
ALIAS(A1,A2);
ALIAS(D1,D2);
ALIAS(D1,D3);

* define the "global" set of scenarios
SET SCEN(A1,A2,D1,D2,D3);
SCEN(A1,A2,D1,D2,D3)=YES;

* define the generator availability fractions
PARAMETER FRAC(I,A1,A2,D1,D2,D3);
FRAC('I1','F1',A2,D1,D2,D3)=1.0;
FRAC('I1','F2',A2,D1,D2,D3)=0.95;
FRAC('I1','F3',A2,D1,D2,D3)=0.3;
FRAC('I2',A1,'F1',D1,D2,D3)=1.0;
FRAC('I2',A1,'F2',D1,D2,D3)=0.8;
FRAC('I2',A1,'F3',D1,D2,D3)=0.0;
* define the demand scenarios

PARAMETER DEMAND(J,A1,A2,D1,D2,D3);
DEMAND('J1',A1,A2,'L1',D2,D3)=900;
DEMAND('J1',A1,A2,'L2',D2,D3)=1000;
DEMAND('J1',A1,A2,'L3',D2,D3)=1300;
DEMAND('J2',A1,A2,D1,'L1',D3)=900;
DEMAND('J2',A1,A2,D1,'L2',D3)=1000;
DEMAND('J2',A1,A2,D1,'L3',D3)=1250;
DEMAND('J3',A1,A2,D1,D2,'L1')=900;
DEMAND('J3',A1,A2,D1,D2,'L2')=1100;
DEMAND('J3',A1,A2,D1,D2,'L3')=1400;

* define marginal probability mass functions
PARAMETER PA1(A1)
	  /F1 0.9
	   F2 0.05
	   F3 0.05 /;
PARAMETER PA2(A2)
	  /F1 0.85
	   F2 0.10
	   F3 0.05 /;
PARAMETER PD(D1)
	  / L1 0.35
	    L2 0.55
	    L3 0.10 /;
* define joint probability mass function (assuming independence)
PARAMETER PROB(A1,A2,D1,D2,D3);

PROB(A1,A2,D1,D2,D3)=PA1(A1)*PA2(A2)*PD(D1)*PD(D2)*PD(D3);

* capacity unit investment cost
PARAMETER KINVEST(I)
	  /I1 400.0
	   I2 350.0 /;

* transportation cost
TABLE CTRANS(I,J)
    J1 J2 J3
I1 4.3 2.0 0.5
I2 7.7 3.0 1.0 ;

* subcontracting cost
SCALAR SUBCON /6000.0/;

* maximum investment
SCALAR BMAXINV /10000.0/;
VARIABLES 
	  Z;
POSITIVE VARIABLES X,Y,S;

EQUATIONS
	UTILITY
	INVCAP
	SUPPLY(I,A1,A2,D1,D2,D3)
	LOAD(J,A1,A2,D1,D2,D3);

* objective is sum of first stage investment cost and expected
* second stage transportation and subcontracting cost
UTILITY .. Z =E= SUM(I,KINVEST(I)*X(I)) + SUM(SCEN,PROB(SCEN)*
		   (SUM((I,J),CTRANS(I,J)*Y(I,J,SCEN))
			   + SUM(J,SUBCON*S(J,SCEN)))) ;

* we can invest no more than BMAXINV units of capacity
INVCAP .. SUM(I, X(I)) =L= BMAXINV ;

* the second stage generation bounds
SUPPLY(I,SCEN) .. SUM(J,Y(I,J,SCEN)) - FRAC(I,SCEN)*X(I) =L= 0.0;

* the second stage demand constraints
LOAD(J,SCEN) .. SUM(I,Y(I,J,SCEN)) + S(J,SCEN) =G= DEMAND(J,SCEN) ;

MODEL CAPACITY /UTILITY, INVCAP, SUPPLY, LOAD/;

SOLVE CAPACITY USING LP MINIMIZING Z;

DISPLAY X.L, Z.L;


SET
	iter 'max benders iterations' /iter1 * iter50/
	dyniter(iter)
;

SCALAR
	lowerbound /-INF/
	upperbound /INF/
;



* Form the master problem
VARIABLES
	ZMASTER {fisrst stage cost}
	theta
	X(I)    {first stage capacity allocation}
	
;

EQUATIONS
	UTILITY
	INVCAP
	OPTCUTS(DYNITER)
;

UTILITY .. ZMASTER =E= SUM(I,KINVEST(I)*X(I)) + SUM(SCEN,PROB(SCEN)*
		   (SUM((I,J),CTRANS(I,J)*Y(I,J,SCEN))
			   + SUM(J,SUBCON*S(J,SCEN)))) ;
* we can invest no more than BMAXINV units of capacity
INVCAP .. SUM(I, X(I)) =L= BMAXINV ;

	
MODEL MASTERPROBLEM / /;

* Form the subproblem
VARIABLES
	ZSUB   {second stage cost}
	Y(I,J,A1,A2,D1,D2,D3) {second stage transportation vars}
	S(J,A1,A2,D1,D2,D3) {second stage subcontracting vars}
	
;

MODEL SUBPROBLEM / /;


* Outer loop through benders iterations
loop (iter$(not done),
	iteration = ord(iter);

*	Inner loop through all subproblems
	loop(s,
		
	)
)



* count is an indicator variable taking on value one if we failed
* to meet demand in a particular scenario and zero otherwise
PARAMETER COUNT(A1,A2,D1,D2,D3);
COUNT(SCEN)=1$(SUM(J,S.L(J,SCEN)) GT 0.0000001);
* intcount tells the probability that we fail to meet demand

SCALAR INTCOUNT;

INTCOUNT=SUM(SCEN,PROB(SCEN)*COUNT(SCEN));

DISPLAY "probability that demand not satisfied = ", INTCOUNT;